#############################################
#Robustness check: Logistic specification
#############################################

rm(list = ls())

#Functions
source("functions.R")

#Packages
installPackageNotFound("stargazer")
installPackageNotFound("survival")
installPackageNotFound("lmtest")

#############################################
#Read in Voter data
#############################################

#All voters
all.voters <- read.csv("voters.csv", stringsAsFactors = F)

#Drop voter pairs in same household as volunteer
all.voters <- all.voters[which(all.voters$same_household == 0), ]

#############################################
#Main regression results
#############################################
#Covariates
covariates <- c("Black", "Hispanic", "Asian", "Male", "Under30",
                "VotePropensity", "Partisanship")

#Data groups
data.list <- list("in_district" = which(all.voters$Volunteer_In_District == 1),
                  "out_district" = which(all.voters$Volunteer_In_District == 0))

#Specifications
itt.base <- paste0("BinaryTurnout~Treatment + strata(pair_id)")
itt.controls <- paste0(itt.base, "+",
                       paste0(covariates, collapse = "+"))
#List of formulas
formula.list <- list("base" = itt.base, "controls" = itt.controls)

#Empty list to hold model results
model.list <- list()
se.list <- list()

for(k in 1:length(data.list)){
  for(i in 1:length(formula.list)){
    #Subset to volunteer group
    this.group <- all.voters[data.list[[k]], ]

    #Treatment effect at the pair level
    this.model <- clogit(formula(formula.list[[i]]), data = this.group,
                         method = "exact")

    #Robust SEs
    this.se <- coeftest(this.model)

    #Store results
    model.list[[names(data.list)[k]]][[i]] <- this.model
    se.list[[names(data.list)[k]]][[i]] <- this.se
  }
}

#############################################
#ITT Latex table
#############################################
title <- "Intention-to-treat Models"

#Full sample model
latex_output <- stargazer(model.list[["in_district"]][[1]],
                          model.list[["in_district"]][[2]],
                          model.list[["out_district"]][[1]],
                          model.list[["out_district"]][[2]],
                          coef = list(coef(model.list[["in_district"]][[1]]),
                                      coef(model.list[["in_district"]][[2]]),
                                      coef(model.list[["out_district"]][[1]]),
                                      coef(model.list[["out_district"]][[2]])
                          ),
                          se = list(se.list[["in_district"]][[1]][, "Std. Error"],
                                    se.list[["in_district"]][[2]][, "Std. Error"],
                                    se.list[["out_district"]][[1]][, "Std. Error"],
                                    se.list[["out_district"]][[2]][, "Std. Error"]
                          ),
                          omit = covariates,
                          covariate.labels = c("Treatment"),
                          omit.stat = "all",
                          add.lines = list(c("Voter Covariates",
                                             "\\multicolumn{1}{c}{N}",
                                             "\\multicolumn{1}{c}{Y}",
                                             "\\multicolumn{1}{c}{N}",
                                             "\\multicolumn{1}{c}{Y}")),
                          column.separate = c(1,1,1,1),
                          dep.var.labels = "",
                          dep.var.caption = "",
                          notes.append = FALSE,
                          column.sep.width = "0pt",
                          no.space = TRUE,
                          title = title,
                          star.cutoffs = c(0.05, 0.01, 0.001),
                          type = "text")

#############################################
#Conditional on assignment order
#############################################
#Empty list to hold model results
model.list <- list()

for(k in 1:length(data.list)){
  for(i in 1:length(formula.list)){
    for(l in 1:4){

      #Subset to volunteer group
      this.group <- all.voters[data.list[[k]], ]
      this.group <- this.group[this.group$assignment_factor == l, ]

      #Treatment effect at the pair level
      this.model <- clogit(formula(formula.list[[i]]), data = this.group,
                           method = "exact")

      #Store results
      model.list[[names(data.list)[k]]][[names(formula.list)[i]]][[l]] <- this.model
    }
  }
}

#Gather coefficients
coef.list <- list(); se.list <- list()
k = 1
for(this.data in c("in_district", "out_district")){
  for(this.formula in c("base", "controls")){
    this.coef <- sapply(1:4,
      function(x)coef(model.list[[this.data]][[this.formula]][[x]])[["Treatment"]])
    this.se <- sapply(1:4, function(x){
      sqrt(vcov(model.list[[this.data]][[
        this.formula]][[x]])[["Treatment", "Treatment"]])
    })
    coef.list[[k]] <- this.coef
    se.list[[k]] <- this.se
    k <- k + 1
  }
}

#############################################
#CATE Latex table
#############################################
title <- "CATE Models"

#Full sample model
latex_output <- stargazer(model.list[["in_district"]][[2]][[1]],
                          model.list[["in_district"]][[2]][[1]],
                          model.list[["out_district"]][[2]][[1]],
                          model.list[["out_district"]][[2]][[1]],
                          coef = coef.list,
                          se = se.list,
                          omit = covariates[4:7],
                          omit.stat = c("all"),
                          covariate.labels = c("Order 1", "Order 2", "Order 3",
                                               "Order 4+"),
                          add.lines = list(c("Voter Covariates",
                                             "\\multicolumn{1}{c}{N}",
                                             "\\multicolumn{1}{c}{Y}",
                                             "\\multicolumn{1}{c}{N}",
                                             "\\multicolumn{1}{c}{Y}")),
                          column.separate = c(1,1,1,1),
                          dep.var.labels = "",
                          dep.var.caption = "",
                          notes.append = FALSE,
                          column.sep.width = "0pt",
                          no.space = TRUE,
                          title = title,
                          star.cutoffs = c(0.05, 0.01, 0.001),
                          type = "text")
